home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / cholesky.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  2KB  |  88 lines
  1. /* choldecomp - Cholesky decomposition  routines.                      */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlsproto.h"
  11. #else
  12. #include "xlsfun.h"
  13. #endif ANSI
  14.  
  15. /*
  16. choldecomp(a, n)
  17.     RMatrix a;
  18.     int n;
  19. {
  20.   int i, j, k;
  21.   
  22.   for (i = 0; i < n; i++) {
  23.     for (k = 0; k < i; k++)
  24.       a[i][i] -= a[k][i] * a[k][i];
  25.     a[i][i] = (a[i][i] > 0) ? sqrt(a[i][i]) : 0.0;
  26.     for (j = i + 1; j < n; j++) {
  27.       for (k = 0; k < i; k++) a[i][j] -= a[k][i] * a[k][j];
  28.       a[i][j] = (a[i][i] > 0.0) ? a[i][j] / a[i][i] : 0.0;
  29.     }
  30.   }
  31.   for (i = 0; i < n; i++)
  32.     for (j = 0; j < i; j++)
  33.       a[i][j] = 0.0;
  34. }
  35. */
  36. /*  macro in xlsdef.h   JKL
  37. static double Max(a, b)
  38.     double a, b;
  39. {
  40.   return(a > b ? a : b);
  41. }
  42. */
  43. void choldecomp(a, n, maxoffl, maxadd)
  44.     RMatrix a;
  45.     int n;
  46.     double maxoffl, *maxadd;
  47. {
  48.   double minl, minljj, minl2;
  49.   int i, j, k;
  50.   
  51.   minl = pow(macheps(), 0.25) * maxoffl;
  52.   minl2 = 0.0;
  53.   
  54.   if (maxoffl == 0.0) {
  55.     for (i = 0; i < n; i++)
  56.       maxoffl = Max(fabs(a[i][i]), maxoffl);
  57.     maxoffl = sqrt(maxoffl);
  58.     minl2 = sqrt(macheps()) * maxoffl;
  59.   }
  60.   
  61.   *maxadd = 0.0;
  62.   for (j = 0; j < n; j++) {
  63.     for (i = 0; i < j; i++) a[j][j] -= a[j][i] * a[j][i];
  64.     
  65.     minljj = 0.0;
  66.     
  67.     for (i = j + 1; i < n; i++) {
  68.       a[i][j] = a[j][i];
  69.       for (k = 0; k < j; k++) a[i][j] -= a[i][k] * a[j][k];
  70.       minljj = Max(fabs(a[i][j]), minljj);
  71.     }
  72.     
  73.     minljj = Max(minljj / maxoffl, minl);
  74.     
  75.     if (a[j][j] > minljj * minljj) a[j][j] = sqrt(a[j][j]);
  76.     else {
  77.       if (minljj < minl2) minljj = minl2;
  78.       *maxadd = Max(*maxadd, minljj * minljj - a[j][j]);
  79.       a[j][j] = minljj;
  80.     }
  81.     
  82.     for (i = j + 1; i < n; i++) a[i][j] /= a[j][j];
  83.   }
  84.   
  85.   for (i = 0; i < n; i++)
  86.     for (j = i + 1; j < n; j++) a[i][j] = 0.0;
  87. }
  88.